
# CHANGE THIS

setwd("/path/to/directory")
myconfidentialfile <- "/path/to/Study/Data/daF3063e.por"

library(haven)
library(plotly)
library(ggplot2)
library(dplyr)
library(scales)
library(rmapshaper)
library(readxl)
library(stringr)
library(sf)
library(tidyverse)
library(gridExtra)
library(sjPlot)
library(estimatr)
library(coefplot)

###############################################################################

### Figure A1

# shape of the municipalities

finland.adm3.shp <- read_sf(dsn=".", layer="kunta1000k_2019", stringsAsFactors = FALSE)

finland.adm3.shp <- as(finland.adm3.shp, "Spatial")

finland.adm3.shp.df <- fortify(finland.adm3.shp, region = "kunta")

# information about refugee centers

VOK <- as.data.frame(read_excel("data_src_Municodes_map.xlsx"))

VOK$id <- str_pad(VOK$id, width=3, side="left", pad="0")

# Let's unite the two above

finland.adm3.shp.df <- merge(finland.adm3.shp.df, VOK, by = "id")

finland.adm3.shp.df <- arrange(finland.adm3.shp.df, order)

highcolor <- "red"
lowcolor  <- "lightgray"

# let's draw a map
p1 <- ggplot(data = finland.adm3.shp.df, aes(x = long, y = lat, group = group))
p2 <- geom_polygon(aes(fill = Refugees_per_capita))
p3 <- scale_fill_gradient(high = highcolor, low = lowcolor, guide = "colorbar", labels = scales::percent_format(accuracy = 1))
p4 <- coord_fixed(1)
p5 <- guides(fill = guide_colorbar(title = "Asylum seekers per capita"))
p6 <- labs(title = "Proportion of asylum seekers per capita", subtitle = "2015-2016")
p7 <- theme_void()

myplot <- p1+p2+p3+p4+p5+p6+p7

ggsave("Figure_A01.pdf", plot = myplot)

###############################################################################

### Figure A2

Municodesfull<- read_dta("data_src_Municodes.dta")

New_data <- subset(Municodesfull, Treated2015 == 1)

mydensity <- ggplot(New_data, aes(x=`refugees_per_capita`)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white") +
  theme_bw() + labs(y = "Density") + labs(x = "Asylum seekers per capita") +  
  geom_density(alpha=.2, fill="#FF6666") 

ggsave("Figure_A02.pdf", plot = mydensity)

###############################################################################

### Figure A3: votes

# first need to calculate vote share of municipalities as defined today, 
# if they have been united with another municipality in the past
# (how municipality mergers affect vote shares)

mergesheet <- read_excel("data_src_fi_mergersofmuncipalities.xlsx")

muni_election_fi <- read_excel("data_src_Munivotes2000-.xlsx",na=c(".","-")) %>%
  select(Year,Municipality,Total,PS,GREENS, KOK, SDP, KESK, VAS, RKP, KD) %>%
  mutate(
         Total = as.numeric(as.character(Total)),
         PS = as.numeric(as.character(PS)),
         GREENS = as.numeric(as.character(GREENS)),
         KOK = as.numeric(as.character(KOK)),
         SDP = as.numeric(as.character(SDP)),
         KESK = as.numeric(as.character(KESK)),
         VAS = as.numeric(as.character(VAS)),
         RKP = as.numeric(as.character(RKP)),
         KD = as.numeric(as.character(KD))
        )  %>%
  left_join(., mergesheet,by=c("Municipality")) %>%
  group_by(Year,New_name) %>%
  summarise(
            PS = sum(PS,na.rm=TRUE),
            GREENS = sum(GREENS,na.rm=TRUE),
            KOK = sum(KOK,na.rm=TRUE) ,
            SDP = sum(SDP,na.rm=TRUE), 
            KESK = sum(KESK,na.rm=TRUE), 
            VAS = sum(VAS,na.rm=TRUE),
            RKP = sum(RKP,na.rm=TRUE), 
            KD = sum(KD,na.rm=TRUE),
            Total = sum(Total,na.rm=TRUE)
           ) %>%
  mutate(
         psshare = PS/Total,
         greenshare = GREENS/Total, 
         KOKshare = KOK/Total, 
         SDPshare = SDP/Total, 
         KESKshare = KESK/Total, 
         VASshare = VAS/Total, 
         RKPshare = RKP/Total, 
         KDShare = KD/Total
        )  %>%
  drop_na(New_name)

muni_election_fi <- muni_election_fi %>% rename(muni = New_name)
muni_election_fi$Year <- as.numeric(as.character(muni_election_fi$Year))

data_src_Municodes <- read_dta("data_src_Municodes.dta")

votes <- left_join(muni_election_fi, data_src_Municodes, by = "muni")

votes <- votes %>% mutate(treat = if_else(Treated2015 == 1 & Year > 2015, 1, 0))

forPT <- subset(votes, Year <= 2012)
forPTwithoutpretreat <- subset(forPT,Pretreated == 0)

# 1/8

CenterCF <- ggplot(forPTwithoutpretreat, aes(x = Year, y = KESKshare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +
  labs(x = "Year") + labs(y = "Center") + labs(colour = "Treated") + theme_bw() +
  scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))

CenterCF

# 2/8

KOKCF <- ggplot(forPTwithoutpretreat, aes(x = Year, y = KOKshare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)
  ) +  
  labs(x = "Year") + labs(y = "National Coalition") + labs(colour = "Treated") + theme_bw() + 
  scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))

KOKCF

# 3/8

SDPCF <- ggplot(forPTwithoutpretreat, aes(x = Year, y = SDPshare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)
  ) +  
  labs(x = "Year") + labs(y = "Social Democrats") + labs(colour = "Treated") + theme_bw() + 
  scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))

SDPCF

# 4/8

PSCF <- ggplot(forPTwithoutpretreat, aes(x = Year, y = psshare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Finns") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
PSCF

# 5/8

GreenCF <- ggplot(forPT, aes(x = Year, y = greenshare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Green League") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))

GreenCF

# 6/8

VASCF <- ggplot(forPTwithoutpretreat, aes(x = Year, y = VASshare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Left Alliance") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))

VASCF

# 7/8

KDCF <- ggplot(forPTwithoutpretreat, aes(x = Year, y = KDShare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Christian Democrats") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
KDCF

# 8/8

SWCF <- ggplot(forPTwithoutpretreat, aes(x = Year, y = RKPshare, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Swedish People's Party") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))

SWCF

# all

mygrid8 <- grid.arrange(CenterCF, KOKCF, SDPCF, PSCF, GreenCF, VASCF, KDCF, SWCF)

ggsave("Figure_A03.pdf", plot = mygrid8)

###############################################################################

### Figure A4:  Parallel trends for municipality characteristics

parallel2000_15_reduced <- read_dta("data_src_parallel2000-15,reduced.dta")

pretreatdropped <- subset(parallel2000_15_reduced, Pretreated==0)

# 1/6

popcf<- ggplot(pretreatdropped, aes(x = year, y = population, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  
  labs(x = "Year") + labs(y = "Population") + labs(colour = "Treated") + theme_bw() + 
  scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
  
popcf

# 2/6

pochcf <- ggplot(pretreatdropped, aes(x = year, y = populationchange, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Population change") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
pochcf

# 3/6

urbcf <- ggplot(pretreatdropped, aes(x = year, y = urban, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Urban density") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
urbcf

# 4/6

forcf <- ggplot(pretreatdropped, aes(x = year, y = foreign, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Share of foreigners") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
forcf

# 5/6

uncf <- ggplot(pretreatdropped, aes(x = year, y = unemployment, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Unemployment") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
uncf

# 6/6

depcf <- ggplot(pretreatdropped, aes(x = year, y = dependencyratio, na.rm = TRUE, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Dependency ratio") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
depcf

# all

mygrid6 <- grid.arrange(popcf, pochcf, urbcf, forcf, uncf, depcf)

ggsave("Figure_A04.pdf", plot = mygrid6)

###############################################################################

### Figure A5:  Parallel trends for placebo regression about environmental values

X2008_2012_2017_LD_TO_USE <- read_dta("data_src_2008_2012_2017_anonym.dta")
pretreatdropped <- subset(X2008_2012_2017_LD_TO_USE, Pretreated==0 & year!=2017) 

environment_PT<- ggplot(pretreatdropped, aes(x = year, y = env_2, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(0.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)) +  labs(x = "Year") + labs(y = "Environment") + labs(colour = "Treated") + theme_bw() + scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black")) + 
  scale_x_continuous(breaks = c(2008,2012))

environment_PT

ggsave("Figure_A05.pdf", plot = environment_PT)

###############################################################################

# Balance tests

### Figure A6

parallel2000_15_reduced <- read_dta("data_src_parallel2000-15,reduced.dta")
cov2015 <-subset(parallel2000_15_reduced, year==2015)

cov2015preatreateddropped <-subset(cov2015, Pretreated == 0)

modelformula1 <- urban ~ Treated2015 + population + populationchange + swedish + highereduc + unemployment + contributionmargin + foreign
Cov_scaled_data <- lapply(cov2015preatreateddropped[, all.vars(modelformula1)], scale) 

Urban <- lm(urban ~ Treated2015, data=Cov_scaled_data)
Population <- lm(population~ Treated2015, data=Cov_scaled_data)
Pop.Change<- lm(populationchange ~ Treated2015, data=Cov_scaled_data)
Swedish <- lm(swedish ~ Treated2015, data=Cov_scaled_data)
Highereduc <- lm_robust(highereduc ~ Treated2015, data=Cov_scaled_data)
Unemployment <- lm(unemployment ~ Treated2015, data=Cov_scaled_data)
Contr.Margin <- lm(contributionmargin ~ Treated2015, data=Cov_scaled_data)
Foreign <- lm(foreign ~ Treated2015, data=Cov_scaled_data)

myfigA6 <- multiplot(Urban, Population, Pop.Change, Swedish, Highereduc, Unemployment, Contr.Margin, Foreign, intercept = FALSE, coefficients = "Treated2015",
          title = "", xlab = "", ylab = "",
          innerCI = 1, outerCI = 2, lwdInner = 1,
          lwdOuter = 0.5,
          by="Model", horizontal=FALSE, fillColor = "white") + 
  theme_bw()

myfigA6

ggsave("Figure_A06.pdf", plot = myfigA6)

###############################################################################

### Figure A7

cov2015 <-subset(parallel2000_15_reduced, year==2015)

cov2015preatreateddropped <-subset(cov2015, Pretreated==0)

# Urban (a)

Cov2015urban <- subset(cov2015preatreateddropped, urban >= 60 & population >= 15000)

modelformula3 <- urban ~ Treated2015 + population + populationchange + swedish + highereduc + unemployment + contributionmargin + foreign
Cov_scaled_data_urban <- lapply(Cov2015urban[, all.vars(modelformula3)], scale) 

Population <- lm(population ~ Treated2015, data=Cov_scaled_data_urban)
Pop.Change<- lm(populationchange ~ Treated2015, data=Cov_scaled_data_urban)
Swedish <- lm(swedish ~ Treated2015, data=Cov_scaled_data_urban)
Highereduc <- lm(highereduc ~ Treated2015, data=Cov_scaled_data_urban)
Unemployment <- lm(unemployment ~ Treated2015, data=Cov_scaled_data_urban)
Contr.Margin <- lm(contributionmargin ~ Treated2015, data=Cov_scaled_data_urban)
Foreign <- lm(foreign ~ Treated2015, data=Cov_scaled_data_urban)
Urban <- lm(urban ~ Treated2015, data=Cov_scaled_data_urban)

myurban <- multiplot(Urban, Population, Pop.Change, Swedish, Highereduc, Unemployment, Contr.Margin, Foreign, intercept = FALSE, coefficients = "Treated2015",
          title = "", xlab = "", ylab = "",
          innerCI = 1, outerCI = 2, lwdInner = 1,
          lwdOuter = 0.5,
          by="Model", horizontal=FALSE, fillColor = "white") + 
  theme_bw()

myurban

ggsave("Figure_A07a.pdf", plot = myurban)

# Rural (b)

Cov2015rural <- subset(cov2015preatreateddropped, urban < 60 & population <15000)

modelformula2 <- urban ~ Treated2015 + population + populationchange + swedish + highereduc + unemployment + contributionmargin + foreign
Cov_scaled_data_rural <- lapply(Cov2015rural[, all.vars(modelformula2)], scale) 

Population <- lm(population ~ Treated2015, data=Cov_scaled_data_rural)
Pop.Change<- lm(populationchange ~ Treated2015, data=Cov_scaled_data_rural)
Swedish <- lm(swedish ~ Treated2015, data=Cov_scaled_data_rural)
Highereduc <- lm(highereduc ~ Treated2015, data=Cov_scaled_data_rural)
Unemployment <- lm(unemployment ~ Treated2015, data=Cov_scaled_data_rural)
Contr.Margin <- lm(contributionmargin ~ Treated2015, data=Cov_scaled_data_rural)
Foreign <- lm(foreign ~ Treated2015, data=Cov_scaled_data_rural)
Urban <- lm(urban ~ Treated2015, data=Cov_scaled_data_rural)

myrural <- multiplot(Urban, Population, Pop.Change, Swedish, Highereduc, Unemployment, Contr.Margin, Foreign, intercept = FALSE, coefficients = "Treated2015",
          title = "", xlab = "", ylab = "",
          innerCI = 1, outerCI = 2, lwdInner = 1,
          lwdOuter = 0.5,
          by="Model", horizontal=FALSE, fillColor = "white") + 
  theme_bw()

myrural

ggsave("Figure_A07b.pdf", plot = myrural)

###############################################################################

### Figure B1

# A

# Balance between candidates that answer and do not answer VAA
Backgroundinfo2017_extended <- read_dta("data_src_Backgroundinfo2017_withmunis.dta")

sex1 <- lm(filledin ~ Sex, data = Backgroundinfo2017_extended )
language1 <- lm(filledin ~ language, data = Backgroundinfo2017_extended )
treated <- lm(filledin~Treated2015, data=Backgroundinfo2017_extended )
age1 <- lm(filledin ~age, data = Backgroundinfo2017_extended )

urban <- lm(filledin ~ Urban_2017, data= Backgroundinfo2017_extended) #for this need to add municovariates2017 to backgroundinfo

##unite all coefficients in one graph##
myplota <- plot_models(treated, age1, sex1, language1, urban, std.est = "std", 
  m.labels = c("Treated", "Age", "Sex", "Mother tongue", "Urban"), 
  auto.label = FALSE, axis.title = "Balance", legend.title = "", colors = "gs")

myplota

ggsave("Figure_B01a.pdf", plot = myplota)

# B

# Balance between candidates in treated and non-treated areas
First_difference_dataset <- read_dta("data_YLE_wide.dta")

First_difference_dataset_pre_treated_dropped <- subset(First_difference_dataset, Pretreated==0)

answeredtwice <- lm(answeredtwice~ Treated2015, data= First_difference_dataset_pre_treated_dropped)

new <- lm(new~Treated2015, data = First_difference_dataset_pre_treated_dropped)

##Add further characteristics from the candidate background info dataset###

Backgroundinfo2017 <- read_dta("data_src_Backgroundinfo2017_withmunis.dta")

Backgroundinfo2017_pre_treated_dropped<- subset(Backgroundinfo2017, Pretreated == 0)

sex <- lm(Sex~Treated2015, data=Backgroundinfo2017_pre_treated_dropped)
language <-lm(language~Treated2015, data=Backgroundinfo2017_pre_treated_dropped)
income <- lm(Income~Treated2015, data=Backgroundinfo2017_pre_treated_dropped)
education<- lm(Education~Treated2015, data=Backgroundinfo2017_pre_treated_dropped)
age <- lm(age ~ Treated2015, data= Backgroundinfo2017_pre_treated_dropped)
##unite all coefficients in one graph##

myplotb <- plot_models(answeredtwice, age, new, sex, language,income,education, 
  std.est = "std", m.labels = c("Answered twice", "Age", "New", "Sex", "Mother tongue", "Income", "Education"), 
  auto.label = FALSE, axis.title = "Balance", legend.title = "", colors = "gs")

myplotb

ggsave("Figure_B01b.pdf", plot = myplotb)

###############################################################################

### Figure E2

# a

muni_economics <- read_dta("data_src_muni_economics.dta")

pretreatdroppedecon <- subset(muni_economics, Pretreated==0)

econpr <- ggplot(pretreatdroppedecon, aes(x = year, y = contribution_margin, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)
              ) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)
              ) +  
  labs(x = "Year") + labs(y = "Contribution margin, euros") + labs(colour = "Treated") +  theme_bw() + 
  scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))
  
econpr

ggsave("Figure_E02a.pdf", plot = econpr)

# b

Tax <- ggplot(pretreatdroppedecon, aes(x = year, y = Taxincome, color = factor(Treated2015))) + 
  stat_summary(fun = mean,
               fun.min = function(x) mean(x) - 1.96*sd(x)/sqrt(length(x)), 
               fun.max = function(x) mean(x) + 1.96*sd(x)/sqrt(length(x)), 
               geom = "pointrange",
               position = position_dodge(.5)) +
  stat_summary(fun = mean,
               geom = "line",
               position = position_dodge(.5)
              ) +  labs(x = "Year") + labs(y = "Taxincome, euros") + labs(colour = "Treated") +  
  theme_bw() + labs(colour = "Treated") +  theme_bw() + 
  scale_color_manual(labels = c("Non-Treated", "Treated"), values = c("red", "black"))

Tax

ggsave("Figure_E02b.pdf", plot = Tax)

###############################################################################

### NOTE: TABLES A3 and A4 don't rely on my own data analysis but the values are taken directly from the file "cbF3063e.pdf"

daF3063e <- read_por(myconfidentialfile)

rural <- subset(daF3063e, BV3==8)
urban <- subset(daF3063e, BV3!=8)

###############################################################################

### Table A5

Q1tab <- table(rural$Q1)
prop.table(Q1tab)

# 1          2          3 
# 0.58720930 0.36627907 0.04651163

Q1taburban <- table(urban$Q1)
prop.table(Q1taburban)

# 1          2          3 
# 0.65786315 0.28091236 0.06122449 

t.test(rural$Q1, urban$Q1) # not significant

Q2tab <- table(rural$Q2)
prop.table(Q2tab)

# 1          2          3 
# 0.20930233 0.69767442 0.09302326 

Q2taburban <- table(urban$Q2)

prop.table(Q2taburban)

# 1          2          3 
# 0.16446579 0.78391357 0.05162065 

t.test(rural$Q2, urban$Q2) # no difference

###############################################################################

### Table A6:  Differences in attitudes between rural and urban people##

# Finland should do more

Q3tab <- table(rural$Q3_1)
Q3taburban<-table(urban$Q3_1)

t.test(rural$Q3_1, urban$Q3_1) #sig

prop.table(Q3tab)
# 1          2          3          4          5 
# 0.16279070 0.36046512 0.20930233 0.17441860 0.09302326 
prop.table(Q3taburban)
# 0.17767107 0.40216086 0.23529412 0.13685474 0.04801921 

# Not in need of protection

Q3tab_2 <- table(rural$Q3_2)
Q3taburban_2 <- table(urban$Q3_2)

t.test(rural$Q3_2, urban$Q3_2) # sig

prop.table(Q3tab_2)

# 1          2          3          4          5 
# 0.19186047 0.30813953 0.30813953 0.13953488 0.05232558 

prop.table(Q3taburban_2)
# 0.14885954 0.26410564 0.33373349 0.17647059 0.07683073 

# Benefits too generous

Q3tab_3 <- table(rural$Q3_3)
Q3taburban_3 <- table(urban$Q3_3)

t.test(rural$Q3_3, urban$Q3_3) # sig

prop.table(Q3tab_3)
# 1          2          3          4          5 
# 0.40697674 0.32558140 0.16860465 0.05232558 0.04651163 

prop.table(Q3taburban_3)

# 0.34093637 0.28571429 0.17647059 0.10204082 0.09483794 

# Young men problematic

Q4tab_1 <- table(rural$Q4_1)
Q4taburban_1 <- table(urban$Q4_1)

t.test(rural$Q4_1, urban$Q4_1)

prop.table(Q4tab_1)

# 1          2          3          4          5 
# 0.30813953 0.29651163 0.19186047 0.11046512 0.09302326 

prop.table(Q4taburban_1)

# 0.25450180 0.31932773 0.19927971 0.13205282 0.09483794 

# Muslims problematic

Q4tab_2 <- table(rural$Q4_2)
Q4taburban_2 <- table(urban$Q4_2)

t.test(rural$Q4_2, urban$Q4_2)

prop.table(Q4tab_2)

# 1         2         3         4         5 
# 0.1686047 0.2616279 0.3197674 0.0872093 0.1627907 

prop.table(Q4taburban_2)

# 0.1272509 0.2701080 0.3265306 0.1608643 0.1152461 

# Never taxes problematic

Q4tab_7 <- table(rural$Q4_7)
Q4taburban_7 <- table(urban$Q4_7)
prop.table(Q4tab_7)

# 1          2          3          4          5 
# 0.39534884 0.37209302 0.05232558 0.13372093 0.04651163 

prop.table(Q4taburban_7)

# 0.29411765 0.34813926 0.07322929 0.21608643 0.06842737 

t.test(rural$Q4_7, urban$Q4_7) #sig

# More crime problematic

Q4tab_3 <- table(rural$Q4_3)
Q4taburban_3 <- table(urban$Q4_3)
prop.table(Q4tab_3)

# 1          2          3          4          5 
# 0.33139535 0.44186047 0.08139535 0.09302326 0.05232558 

prop.table(Q4taburban_3)
# 0.25090036 0.46578631 0.09603842 0.14405762 0.04321729 
t.test(rural$Q4_3, urban$Q4_3)

# Dark skin problematic

Q4tab_4 <- table(rural$Q4_4)
Q4taburban_4 <- table(urban$Q4_4)
prop.table(Q4tab_4)
# 1          2          3          4          5
# 0.06976744 0.13953488 0.63372093 0.09302326 0.06395349 
prop.table(Q4taburban_4)
# 0.05162065 0.10444178 0.68787515 0.10564226 0.05042017 
t.test(rural$Q4_4, urban$Q4_4)

# Too well to do problematic

Q4tab_5 <- table(rural$Q4_5)
Q4taburban_5<-table(urban$Q4_5)
prop.table(Q4tab_5)

# 1         2         3         4         5 
# 0.1627907 0.3197674 0.2267442 0.1802326 0.1104651 
prop.table(Q4taburban_5)
# 0.11404562 0.26170468 0.29051621 0.25570228 0.07803121 

t.test(rural$Q4_5, urban$Q4_5)

# No integration problematic 

Q4tab_6 <- table(rural$Q4_6)
Q4taburban_6<-table(urban$Q4_6)
prop.table(Q4tab_6)
# 1          2          3          4          5 
# 0.22674419 0.53488372 0.09302326 0.08720930 0.05813953 
prop.table(Q4taburban_6)
# 0.21368547 0.47058824 0.12725090 0.15246098 0.03601441
t.test(rural$Q4_6, urban$Q4_6)

